% script to plot dependence of file size, computation time & peak
% wavelength on mesh cell size
figurepath='X:\josh-phd-code\cpgs\cpgs\4_ResultsAndDiscussion\ResultsAndDiscussionFigs\';
pathname='Y:\2010-11-10-bowties-with-symmetric-boundaries\';
basename='73_09_60_30_16_m';
meshsizes=[0.25:0.05:1];
fontsize=10;
simtimes=zeros(size(meshsizes));
filesizes=[15090381185,... % 0.25
           10507119241,... % 0.3
           7791218665,...  % 0.35
           5983616569,...  % 0.4
           4720726681,...  % 0.45
           3832092121,...  % 0.5
           3204415449,...  % 0.55
           2692143161,...  % 0.6
           2312893489,...  % 0.65
           2001902689,...  % 0.7
           1752410849,...  % 0.75
           1546901614,...  % 0.8
           1365553473,...  % 0.85
           1221631742,...  % 0.9
           1102706022,...  % 0.95
           992223369];     % 1
%           2479841334,...  % 0.625           
%           2168021534,...  % 0.675
%           1864099078,...  % 0.725
%           1652532214,...  % 0.775
%           1458934094,...  % 0.825
%           1284698182,...  % 0.875
%           1177325798,...  % 0.925
%           1049266710,...  %0.975

for ii=1:length(meshsizes)
    f=fopen([pathname basename num2str(100*meshsizes(ii)) '_p0.log']);
    while ~feof(f)
        p='time to run.*: (\d+)';
        m=regexp(fgetl(f),p,'tokens');
        if ~isempty(m)
            simtimes(ii)=simtimes(ii)+str2double(m{1});
        end
    end
    fclose(f);
end

mask=logical((abs(meshsizes-0.4)<0.04) +...
             (abs(meshsizes-0.6)<0.04) +...
             (abs(meshsizes-0.85)<0.04) +...
             (abs(meshsizes-0.25)<0.04));
simtimes(mask)=simtimes(mask)/2;


for ii=1:length(meshsizes)
    n=[pathname basename num2str(100*meshsizes(ii)) '_data.txt'];
    data=load(n);
    lambda=data(:,1);
    sigma_scat=data(:,end-2);
    pk=lambda(sigma_scat==max(sigma_scat));
    l_int=(pk-15):0.0001:(pk+15);
    s_int=spline(lambda,sigma_scat,l_int);
    p2(ii)=l_int(s_int==max(s_int));
end

%%

figure(1)
clf(1)
hold on
ax(1) = subplot(1,2,1);
set(ax(1),'xtick',[],'ytick',[],'box','on')
ax(2) = axes('Units',get(ax(1),'Units'), ...
    'Position',get(ax(1),'Position'),'Parent',get(ax(1),'Parent'));
ax(3) = axes('Units',get(ax(2),'Units'), ...
    'Position',get(ax(2),'Position'),'Parent',get(ax(2),'Parent'));

%[ax h1 h2]=plotyy(h,meshsizes,simtimes/3600,meshsizes,filesizes(filesizes~=0)/(1024^3));
h1=plot(ax(2),meshsizes,simtimes/3600);
h2=plot(ax(3),meshsizes,filesizes(filesizes~=0)/(1024^3));
hold off
set(ax(3),'YAxisLocation','right','Color','none', ...
    'XGrid','off','YGrid','off','Box','off', ...
    'HitTest','off');
set(h1,'color','b')
set(h2,'color','r')
set(h1,'Marker','x')
set(h2,'Marker','+')
set(ax(2),'YColor','k')
set(ax(3),'YColor','k')
set(get(ax(2),'YLabel'),'String','simulation time /hrs')
set(get(ax(3),'YLabel'),'String','file size /Gb')
ax=[ax(1) ax(2) ax(3)];
set([ax get(ax(2),'Xlabel') get(ax(2),'Ylabel') get(ax(3),'Ylabel')],...
    'Fontsize',fontsize,'FontName','Times New Roman')
set(get(ax(2),'xlabel'),'string','mesh size /nm');
legend([h1 h2],'time','file size')
set(ax(2),'YTick',(0:4:round(max(simtimes/3600))),...
          'YTickLabel',num2str((0:4:round(max(simtimes/3600)))'));
set(ax(3),'YTick',0:2:max(filesizes/(1024^3)),...
          'YTickLabel',num2str((0:2:max(filesizes/(1024^3)))'));
set(ax(3),'xtick',[])
set(ax(2:3),'box','off','color','none')
%ax(4)=axes('color','none','box','on','xtick',[],'ytick',[]);
%ax=[ax ax(4)];
set(ax,'position',get(ax(1),'position'))
%%
axnum=4;
figure(1)
ax(axnum)=subplot(1,2,2);
plot(ax(axnum),meshsizes,p2,'g.-')
xlabel('mesh size /nm')
ylabel('\lambda_{max} /nm')
hold on
%plot(0.25:0.001:1,spline(meshsizes,p2,0.25:0.001:1),'r')
hold off
set([ax(axnum) get(ax(axnum),'Xlabel') get(ax(axnum),'Ylabel')],...
    'Fontsize',fontsize,'FontName','Times New Roman')
%%
set(ax,'units','normalized')
set(ax(1:3),'position',[0.07 0.15 0.37 0.82])
set(ax(4),'position',[0.605 0.15 0.38 0.82])


%%
set(gcf,'PaperSize',[16,6],'PaperPositionMode','auto',...
        'units','centimeters','paperposition',[0 0 16 6],'position',[0 2 16 6]*0.8);
%set(gca,'units','normalized');
%set(gca,'position',[0.1 0.075 0.85 0.85]);
print('-dpdf',[figurepath 'FDTDmeshsizes']);
